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ABSTRACT 

We perform a comparative numerical hydrodynamics study of embedded protostellar disks formed 
as a result of the gravitational collapse of cloud cores of distinct mass (M c i=0.2-1.7 M©) and ratio 
of rotational to gravitational energy (/?=0. 0028-0. 023). An increase in M c \ and/or f3 leads to the for- 
mation of protostellar disks that are more susceptible to gravitational instability. Disk fragmentation 
occurs in most models but its effect is often limited to the very early stage, with the fragments being 
either dispersed or driven onto the forming star during tens of orbital periods. Only cloud cores with 
high enough M c \ or /? may eventually form wide-separation binary /multiple systems with low mass 
ratios and brown dwarf or sub-solar mass companions. It is feasible that such systems may eventually 
break up, giving birth to rogue brown dwarfs. Protostellar disks of equal age formed from cloud cores 
of greater mass (but equal f3) are generally denser, hotter, larger, and more massive. On the other 
hand, protostellar disks formed from cloud cores of higher f3 (but equal M c \) are generally thinner 
and colder but larger and more massive. In all models, the difference between the irradiation temper- 
ature and midplane temperature AT is small, except for the innermost regions of young disks, dense 
fragments, and disk's outer edge where AT is negative and may reach a factor of two or even more. 
Gravitationally unstable, embedded disks show radial pulsations, the amplitude of which increases 
along the line of increasing M c \ and j3 but tends to diminish as the envelope clears. We find that 
single stars with a disk-to-star mass ratio of order unity can be formed only from high-/3 cloud cores, 
but such massive disks are unstable and quickly fragment into binary /multiple systems. A substan- 
tial fraction of an embedded disk, especially its inner regions, spiral arms and dense clumps, may be 
optically thick, leading potentially to observational underestimates of disk masses in the embedded 
phase of star formation. 



Subject headings: circumstellar matter — planetary systems: protoplanetary disks 
— ISM: clouds — stars: formation 
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1. INTRODUCTION 

The early phases of the evolution of a protostellar disk, 
starting from its formation and ending with the clear- 
ing of a natal cloud core, determine the course along 
which a young stellar object (YSO) will evolve later in 
the T Tauri phase. Vigourous gravitational instabil- 
ity that develops in protostellar disks during the em- 
bedded phase of star formation (EPSF) serves to limit 
disk masses, effectively setting an upper limit on the 
mass accretion rates in the T Tauri phase and flatten- 
ing the ac cretion rate-stellar mass relation for solar- 
type stars ([Vorobvov fc Basul l2009at ). Disk fragmenta- 
tion triggered by gas infall onto the disk from the col- 
lapsing cloud core could lead to the formation of gi- 
ant planets or brown dwarfs, thus shaping the subse- 
quent phy sical properties of stellar and planetary sys- 
tems (e.g. lBolevll2QQ9t iBolev et al.l 120091: iMachida et all 
120101 : IVorobvov fe Basul l2010at ). The large spread in 
mass accretion rates inferred for embedded YSOs in 
Perseus, Serpens, a nd Ophiuchus star-forming regions by 
lEnoch et all (2009) can be accounted for by gravitational 
instabilit y and fragme ntation of protostellar disks in the 
EPSF (|Vorobvovll2009bl ) . In ad dition, episodic accretion 
caused by disk fragmentation ([ Vorobvov fc Basul 12005. 
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I2006L l2010bh can explain the observed luminosity spread 
i n H-R diagrams of several-Myr-old star- forming regions 
([Baraife et al.l I2009D and resolve the lon g standing "lu- 
minosity problem" ([Dunham et al.ll2010h . whereby pro- 
tostars are underluminous compared to the accretion lu- 
minosity expected from analytical collapse calculations. 

In spite of the pivotal role of the EPSF, our knowl- 
edge of protostellar disks in this phase is embarrass- 
ingly inadequate owing to difficulties with both obser- 
vations and modeling. Observing this phase has been 
difficult because disks are only visible at wavelengths 
where resolution is poor. The added complication of 
the envelope structure in the embedded systems makes 
an extraction of disk parameters from interferometer 
data extraordinarily difficult with modern facilities. No 
wonder that most observations of disk properties have 
been focused on T Tauri disks (e.g. [A keson et al. 2002; 



Kita mura et al-ll2QQ2t IVicente fc A lves 200 51: iPietuet al.l 
20071: lAn drews fc Williams! 120071 : lAndrews et al.l I200S 



Isella et all I2009L and many others) and considerably 



fewer attempted studies of the embedded di s ks have been 
done so far (e.g. lAndrews fc Williamsl 120051 : lEisner et al.l 
l2QQ5Uj0rgensen et al.ll2009f ). 

Self-consistent numerical simulations of protostellar 
disks in the EPSF is no less difficult than observations 
due to vastly changing spatial and temporal scales in- 
volved. The matter is that it is not sufficient to just 
consider an isolated system with some presumed disk-to- 
star mass ratio. A self-consistent treatment of the inter- 
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action of the star/disk system with the natal cloud core 
is of considerable importance for the disk physics and 
this inevitably requires solving for a much larger spa- 
tial volume than in isolated systems. As a consequence, 
most information about the early evolution of embed- 
ded protostellar disks was drawn from one-dimensional 
numerical studies of thin axisymmteric disks with some 
approximate prescription for disk gravity and, in some 
cases, with a solution for the vertical dis k structure in 
the so-called 1+1 D approximation (e.g. iLin fc Pringld 
19901: iBell fc Lid 11994 iNakamo tofc Nakagawal 119941: 
Hueso fc Guillotl 12005b IVisser et al.l 12009b iRice at al.l 
2Q1Q[ ). Alternatively, analytic models of protostellar 
disks subject to intense mass loading have been con- 
structed, which cover a wide parameter space and can 
yield important information than can later be com- 
pared against more focused numerical hydrodynam- 
ics simulations (e.g .IKratter et all 120081 : IClarkd 120091 : 
iKrumholz fc Burkertll2010l )~ 

Some interesting insight into the early disk evolution, 
accretion history, and velocity structure of the infalling 
gas was also obtained in two-dime nsional numerical simu- 
lations of axisymmetric disks fe.g.lYorke fc Bodenheimerl 
[T999l : lBoss fc Hartma nn 2001: IZhu et a l. 2009). Unfortu- 
nately, such simulations cannot describe self-consistently 
mass and angular momentum transport due to gravi- 
tational torques, which cannot develop in axisymmetric 
disks. To resolve this problem, some sort of effective vis- 
cosity is often invoked. An alternative two-dimensional 
approach involves solving for the evolution of protostel- 
lar disks in the thin-disk approximation, which allows 
for an accurate treatment of gravitational tor ques (e.g. 
iNelson et all 120001 : Uohnson fc Gammid I2003D . These 
studies revealed the importance of disk cooling for the 
development of gravitational instability and fragmenta- 
tion but suffered from the lack of disk interaction with a 
natal cloud core in the EPSF. 

An imp ortant step forward was done in a se ries of 
papers by IVorobvov fe Basul (|2005l . I2006L l2009al lbh and 
IVorobvov] ( 2009al lbL l2010f ) who employed two-dimensional 
numerical simulations in the thin-disk approximation 
with an accurate treatment of disk self-gravity to model 
self-consistently the formation and evolution of protostel- 
lar disks. These studies have demonstrated that mass 
loading onto the disk from the infalling envelope can 
bring about qualitatively new phenomena such as the 
burst mode of accretion (IVorobvov fc Basul I2005L I2006L 
l201Qbh and formation of giant planets on wide orbits 
(|Vorobvov fc Basul l2010a| ). Such two-dimensional sim- 
ulations allow us to consider the long-term evolution of 
a large number of self-consistently formed protostellar 
disks at the expense of the detailed disk vertical struc- 
ture. 

First fully three-dimensional simulations, though 
suffering from insufficient numerical resolution, have 
revealed important information about the verti- 
cal structure of the disks and their susceptibility 
to gravitational instability and fragmentation (e.g. 
lLaughlin fc Bode nheimer 19941 : iBurkert et al.lll997UBossl 
119981 : iPickett et al.ll2000D . With the advent of more pow- 
erful computers and techniques, new and more detailed 
three-dimensional numerical simulations of the forma- 
tion and evolution of protoste l lar disks have started to 
emerge (|Krumholz et al.l 120071 : iBolevI 120091 : iBolev et al.l 



120091 : IKratter et al.ll2010al iMachida et al.ll201Q[ ). Unfor- 
tunately, due to still enormous computational load in- 
volved, these studies often suffer from restricted param- 
eter space, limited temporal and spatial resolution, sim- 
plified treatment of the gas infall onto the disk, and/or 
the lack of detailed thermal physics. In addition, due to 
severe resolution limitations, 3D numerical simulations 
often resort to the use of sink particles to trace forming 
stellar and planetary objects, which inevitably involves 
introducing additional free parameters into the model. 

This is the second paper in a series addressing the phys- 
ical propert ies of embedded protoste llar disks. In the 
first paper ([Vorobvov fc Basul l201Qbh , we have mainly 
focused on the burst mode of accretion that developes 
in fragmenting disks. In this paper, we perform a com- 
prehensive numerical hydrodynamics study of the forma- 
tion and long-term evolution of protostellar disks formed 
from cloud cores of distinct mass and rotational energy. 
We make use of the thin-disk approximation comple- 
mented by a detailed calculation of the disk thermal bal- 
ance. Several models were scrutinized to derive the de- 
tailed disk structure and time evolution. Our numerical 
studies can serve as a framework for the future high- 
resolution observations of embedded disks using such 
ground-b ased or space-base d facilities as ALMA or Mil- 
limetron (jWild et al.l l2QQ9h . The paper is organized as 
follows. In § [2j we review our numerical model and ini- 
tial setup. The structure and evolution of protostellar 
disks fromed from cloud cores of distinct mass and ro- 
tational energy are presented in § [3] and § [H The main 
conclusions are summarized in § [5j 

2. MODEL DESCRIPTION 

The main concept s of our numerical approa ch are ex- 
plained in detail in IVorobvov fc Basul (|2010b[ ) and, for 
the reader's convenience, are briefly reviewed below. We 
start our numerical integration in the pre-stellar phase, 
which is characterized by a collapsing starless cloud core, 
continue into the embedded phase of star formation, dur- 
ing which a star, disk, and envelope are formed, and ter- 
minate our simulations in the T Tauri phase, when most 
of the envelope has accreted onto the forming star /disk 
system. In the EPSF, the disk occupies the innermost 
region of our numerical grid, while the larger outer part 
of the grid is taken up by the infalling envelope, the latter 
being the remnant of the parent cloud core. This ensures 
that the protostellar disk is not isolated in the EPSF but 
is exposed to intense mass loading form the envelope. In 
addition, the mass accretion rate onto the disk M e nv is 
not a free parameter of the model but is self-consistently 
determined by the gas dynamics in the envelope. 

We introduce a "sink cell" at r sc = 6 AU and impose 
a free inflow inner boundary condition and free outflow 
outer boundary condition so that that the matter is al- 
lowed to flow out of the computational domain but is 
prevented from flowing in. We monitor the gas surface 
density in the sink cell and when its value exceeds a crit- 
ical value for the transition from isothermal to adiabatic 
evolution, we introduce a central point-mass star. In the 
subsequent evolution, 90% of the gas that crosses the in- 
ner boundary is assumed to land onto the central star 
plus the inner axisymmetric disk at r < 6 AU. This in- 
ner disk is dynamically inactive; it contributes only to 
the total gravitational potential and is used to secure a 
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smooth behaviour of the gravity force down to the stellar 
surface. The other 10% of the accreted gas is assumed 
to be carried away with protostellar jets. The latter are 
triggered only after the formation of the central star. 

2.1. Basic equations 

We make use of the thin-disk approximation to com- 
pute the gravitational collapse of rotating, gravitation- 
ally unstable cloud cores. This approximation is an ex- 
cellent means to calculate the evolution for many orbital 
periods and many model parameters. The basic equa- 
tions of mass, momentum, and energy transport are 



as 

dt 

d 



-V p • (Ev p ) , 



(1) 



- (Et7 p ) + [V • (Xv p v p )] p = - V P V + E g p + (2) 

+ (v-n) p , 



d 6 , v 



-V(V p .v p )-A + T + (Vv),:Il p 



(3) 

where subscripts p and p r refers to the planar components 
(r, 0) in polar coordinates, E is the mass surface density, 

e is the internal energy per surface area, V = J*_f z Pdz 
is the vertically integrated form of the gas pressure P, 
Z is the radially and azimuthally varying vertical scale 
height determined in each computational cell using an 
assumption of local hydrostatic equilibrium, v p = v r r + 

v^cj) is the velocity in the disk plane, g p = g r r + g^cj) 
is the gravitational acceleration in the disk plane, and 
V p = rd/dr + (j)r~ 1 d/d(j)is the gradient along the planar 
coordinates of the disk. 

Viscosity enters the basic equations via the viscous 
stress tensor II expressed as 



n = TZv ( Vv - -(V • v)e 



(4) 



where e is the unit tensor. We note that we take no 
simplifying assumptions about the form of II apart from 
those imposed by the adopted thin-disk approximation. 
We parameterize the magnitude of kinematic viscosity v 
using a modified form of the a-prescription 



(5) 



where Cg = 7P/E is the square of effective sound speed 
calculated at each time step from the model's known V 
and E. The function J- a (r) = 2ir~ 1 tan -1 [(jd/r) 10 ] is 
a modification to the usual a-prescription that guaran- 
tees that the turbulent viscosity operates only in the disk 
and quickly reduces to zero beyond the disk radius ra- 
in this paper, we use a spatially and temporally uni- 
form a, wi th its value set to 5 x 1Q ~ 3 based on our re- 
cent work ([Vorobvov fc Ba su 2009b). Our adopted value 
of a are in agreemen t with the mean value inferred by 
I Andrews et all ([2009) for a large sample of protostellar 
disks in the Ophiuchus star-forming region. 

Equation (j3j) for the internal energy (per surface area) 
transport includes compressional heating V (Vp • v p ), ra- 
diative cooling A, heating due to stellar/background ir- 
radiation T, and viscous heating (S7v) pp > : TI PP > . The 



cooling function is determined using the diffusion approx- 
imation of the vertical radiation tra nsport in a one-zone 
mode l of the vertical disk structure ([Johnson fc Gammid 
2003) 



A = F c cf T 4 



1 



(6) 



where a is the Stefan-Boltzmann constant, T is the mid- 
plane temperature of gas, and T Q = 2 + 20 tan -1 (r) / (3ir) 
is a function that secures a correct transition between 
the cooling function (from both surfaces of the disk) in 
the optically thick regime A t hick = 16<rT 4 /3r and the 
optically thin one A t h m = 2aT 4 r. We use frequency- 
integrated opacities of iBell fc Linl ([l994f ) 
The heating function is expressed as 



r = T c aT, 



-4 _ 

irr 2 



(7) 



where T- lvv is the irradiation temperature at the disk sur- 
face determined by the stellar and background black- 
body irradiation as 



PirrW 



(8) 



where Xb g is the uniform background temperature (in our 
model set to the initial temperature of the natal cloud 
core) and F\ TT (r) is the radiation flux (energy per unit 
time per unit surface area) absorbed by the disk surface 
at radial distance r from the central star. The latter 
quantity is calculated as 



Pirr(r) 



Arr^2 COS J h 



(9) 



where 7i rr is the incidence angle of radiation arriving 
at the disk surface at radial distance r and A irr = 
M c \/ (M env + M c \) is a time-dependent factor that ac- 
counts for the attenuation of stellar radiation by the en- 
velope with mass M env in three-dimensional disks and 
takes values from > 0.5 at the disk formation to w 1.0 at 
the end of the EPSF and beyond 3 . The stellar luminos- 
ity L* is the sum of the accretion luminosity a ccr — 
GM*M /2r* arising from the gravitational energy of ac- 
creted gas and the photospheric luminosity L* )P h due to 
gravitational compression and deuterium burning in the 
star interior. The stellar mass M* and accretion rate 
onto the star M are determined self-consistently dur- 
ing numerical simulations via the amount of gas pass- 
ing through the sink cell. The stellar ra dius r* is calcu- 
lated using an approximation formula of iPalla fc Stahlerl 
(|l99lf h modified to t ake into account the formation of the 
first molecular core (Masunaga fc Inutsukall2000f ). The 
photospheric luminosity L*, p h is taken from the pre-main 
sequence tra cks for the low-mass stars and brown dwarfs 
calculated b ylD'Antona fc Mazitel fl (11997ft . More details 
are given m IVorobvov fc Basu! (|2010bf ). 

Viscous heating operates in the disk interior and is 
calculated using the standard expression (Vv) pp / : H pp ' . 

3 Although most of the envelope material should be landing onto 
the disk outer edg e, a smaller f raction may still be falling onto the 
disk inner regions ( Visser et al. 2009) and this material may inter- 
cept some of the stellar radiation in the EPSF. The effect of this 
attenuation factor on the disk evolution is however not significant, 
since its value quickly approaches unity. For instance in model 3, 



at t ■ 



i 0.75 at t ■ 
: 0.3 Myr. 



= 0.1 Myr after the disk formation and A- u 



! 0.95 
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TABLE 1 

Model cloud core parameters 



Model 


M cl 


/3 


^0 


ro 


So 


1 


0.2 


5.6 x 10- 3 


6.2 


445 


0.28 


2 


0.85 


5.6 x 10- 3 


1.5 


1885 


0.066 


3 


1.7 


5.6 x 10- 3 


0.7 


3770 


0.033 


4 


0.85 


2.8 x 10- 3 


1.0 


1885 


0.066 


5 


0.85 


2.3 x 10~ 2 


2.9 


1885 


0.066 



Note. — All masses are in Mq, distances in AU, surface densi- 
ties in g cm -2 , and angular velocities in km s _1 pc _1 . 



Heating due to shock waves is taken into account via 
compressional heating V (V p • v p ) and artificial viscosity, 
the latter impl emented in the code using the standard 
prescription of iRichtmver fc Morton! (|1957| ). The verti- 
cally integrated gas pressure V and internal energy per 
surface area e are related via the ideal equation of state 
V = (7 — 1) e, with the ratio of specific heats 7 = 7/5. 

Equations (PQ)-(|3]) are solved using the method of finite 
differences with a time explicit solution procedure in po- 
lar coordinates (r, 0) on a numerical grid with 512 x 512 
grid zones. The advection is treated using the van Leer 
interpolation scheme. The update of the internal en- 
ergy per surface area e due to cooling A and heating 
T is done implicitly using the Newton-Raphson method 
of root finding, complemented by the bisection method 
where the Newton-Raphson iterations fail to converge. 
The viscous heating and force terms in Equations (|2]) 
and (|3]) are implemented in the code using an explicit 
finite-difference scheme, which is found to be adequate 
for a < 0.01. The radial points are logarithmically 
spaced. The innermost grid point is located at the posi- 
tion of the sink cell r sc = 6 AU, and the size of the first 
adjacent cell varies in the 0.07-0.1 AU range depending 
on the cloud core size. This corresponds to the radial 
resolution of Ar=l. 1-1.6 AU at 100 AU. 



tional energy defined as 

^out 

E rot = 2tt J ra^rdr, (12) 

Tout 

^grav = -27r J rg^rdr, (13) 

respectively, where a c = Vt 2 r is the centrifugal accelera- 
tion. The parameters of these models are listed in Ta- 
ble [U We note that cloud cores with the adopted values 
of M c \ are expected to form solar- and sub-solar type 
sta rs and the values o f /3 lie within the limits inferred 
by ICaselli et al.l (|2002l ) for dense molecular cloud cores, 
P = (10- 4 -0.07). 

In the following text, we consider separately models 
with distinct M c \ but similar j3 (models 1, 2, and 3) 
and distinct j3 but similar M c \ (models 2, 4, and 5). In 
particular, the parameter f3 is varied by changing the 
value of VLq only. These two sets of models have been 
chosen to single out the effects of increasing cloud core 
masses, from one hand, and increasing rotational energies 
of cloud cores, from the other hand. We note that the 
original cloud core mass and rotation rate are expected 
to have a major impact on the evolution of a protostellar 
disk because they determine the centrifugal radius and, 
consequently, the disk mass and stability properties. 

3. DISK STRUCTURE AND EVOLUTION ALONG THE 
LINE OF INCREASING CLOUD CORE MASSES 

We start by considering the properties of protostellar 
disks formed from cloud cores of increasingly higher mass 
M c \. In particular, model 1 is characterized by M c \ = 
0.2 M Q , while models 2 and 3 have M c \ = 0.85 M and 
M c \ = 1.7 Mq, respectively. The ratio f3 is identical for 
these models and is set to 5.6 x 10 -3 . 



2.2. Initial conditions 

Initially isothermal (Ti n i t = Tb g = 10 K) cloud cores 
have surface densities E and angular velocities Q typical 
for a collapsing, a xisymmetric, magnetically supercritical 
core dBasulll997[ ) 




where Qo is the central angular velocity and ro is the 
radius of central near-constant-density plateau defined 
as ro = V^4Cg/(7rGHo), with the initial positive density 
enhancement A set to 1.2 throughout the paper. Model 
cores are characterized by a distinct ratio r ou t/^o = 6 
in order to generate gravitationally unstable truncated 
cores of similar form, where r out is the cloud core's outer 
radius. 

For the in-depth analysis of circumstellar disk struc- 
ture and evolution, we have chosen five typical model 
cloud cores with varying initial cloud core masses M c \ 
and ratios j3 = £?rot/|^grav| of the rotational to gravita- 



3.1. Face-on disk images 

Figure [T] presents gas surface densities (in g cm -2 , log 
units) in model 1 (top row), model 2 (middle row) and 
model 3 (bottom row) at four distinct times after the disk 
formation (from left to right): £=0.1 Myr, £=0.3 Myr, 
£=0.5 Myr, and t = 0.7 Myr. We note that the spatial 
scale is different for model 1 (400 x 400 AU), model 2 
(800 x 800 AU), and model 3 (1100 x 1100 AU). Super- 
imposed on the gas surface density images are the gas 
velocity fields (arrows) and yellow contour lines. The 
latter delineate disk regions that are characterized by 
the frequency- averaged optical depth r > 0.1. Optical 
depth effects are expected to become important in these 
regions. 

It is evident that the disk appearance in Figure Q] 
changes significantly along both the line of increasing 
cloud core mass (from top to bottom) and the line of 
growing disk age (from left to right). The disk in the 
M c \ = 0.2 Mq model 1 (top row) demonstrates a weak, 
diffuse spiral structure at t = 0.1 Myr and becomes 
nearly perfectly axisymmetric in the subsequent evolu- 
tion. On the other hand, the disk in the M c \ = 0.85 Mq 
model 2 (middle row) exhibits a rich, well-defined spiral 
structure at t = 0.1 Myr, which weakens and diffuses 
with time but still can be traced in the 0.5-Myr-old disk. 
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Fig. 1. — Gas surface density (in g cm" , log units) in three models with distinct cloud core mass M c \ and at distinct times t after 
the disk formation. In particular, the top row shows disk images in the M c \ = 0.2 Mq model, while the middle and bottom rows present 
those of the M c \ = 0.85 Mq and M c \ = 1.7 Mq models, respectively. Columns from left to right correspond to t = 0.1 Myr, t = 0.3 Myr, 
t = 0.5 Myr, and t = 0.7 Myr, respectively. Superimposed on the disk images are the gas velocity fields (arrows) and yellow contour lines 
delineating disk regions with the frequency-integrated optical depth r > 0.1. 



Fragmentation of the spiral arms (small red clumps in 
the middle- left panel) occurs at £=0.1 Myr after the disk 
formation. However, no survived fragments are seen at 
later times due to an efficient migration mechanism 4 that 
operates in the EPSF and drives the fragments into the 
inner disk regions and through the sink cell, causing ac- 
cretion and luminosity bursts ([Vorobvov fc Basul |2005, 
2006. l2010bh . An apparent lopsidedness is evident in the 
0.3-Myr-old and 0.5-Myr-old disks but seems to diminish 
at later times. The gas velocity field indicates consider- 
able non-circular motions within the disk, suggesting the 
existence of local deviations from a Keplerian rotation. 
In addition, notable non-axisymmetric variations in the 
disk shape suggest that an axisymmetric description of a 
circumstellar disk in the early evolution stage is unjusti- 
fied, at least on time scales of the order of 0.5 Myr after 
the disk formation. 

The M c \ = 1.7 model 3 (bottom row in Figure [1]) is 
most extreme and demonstrates vigourous gravitational 
instability and fragmentation. Two well-defined frag- 
ments possessing counterrotating minidisks and a few 
smaller ones are seen in the 0.3-Myr-old disk (note that 

4 See animation of this quick migration at 
www.astro.uwo.ca/~vorobyov (animations: burst mode of 
accretion). 



the disk rotates counterclockwise). The counterrotation 
can clearly be seen in Figure [2j which shows the resid- 
ual velocity field superimposed on the gas surface den- 
sity in model 3 at t = 0.3 Myr. The residuals are ob- 
tained by subtracting the Keplerian rotation from the 
total gas velocity, thus accentuating the gas motion in a 
local frame of reference moving with the disk. Consider- 
able non-Keplerian motions are also evident near/at the 
spiral arms. 

Such counterotating structures may owe their existence 
to differential rotation of fragmenting spiral arms. When 
a fragment forms near corotation 5 , it captures gravita- 
tionally some of the neighboring material, which receives 
a counterrotating twist around the forming fragment due 
to the fact that the inner parts of the arm rotate around 
the central star faster than the outer ones. It is impor- 
tant to note that this mechanism is not univers al and, 
as found in other runs (jVorobvov fc Basu 2010b), some 
fragments may also possess corotating minidisks. It is 
feasible that the direction of rotation is determined by 
local conditions at the position of the forming fragment 
(e.g. spiral arm structure and density, residual veloc- 
ity field, etc.), in which case it is difficult to predict the 

5 Fragmentation of spiral arms occurs preferentially near coro- 
tation where distortion due to differential rotation is minimal. 
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Fig. 2. — Residual velocity field (total velocity field minus Keple- 
rian rotation) superimposed on the gas surface density in model 3 
at t=0.3 Myr after the disk formation. The horizontal arrow in 
the right-upper corner has a dimension of 5 km s — 1 . In the sim- 
ulation, the disk rotates counterclockwise and two fragments with 
disk rotating clockwise are clearly visible. 

outcome. We plan to study this phenomenon in greater 
detail in a future paper. 

In the 0.5-Myr-old disk, most of the fragments have 
either been destroyed or absorbed by the central star. 
However, one fragment has survived and a binary system 
with the smaller companion possessing a sizeable coun- 
terotating disk has emerged. The separation between 
the companions is about 1000 AU and the masses of the 
primary /secondary are approximately M* = 0.95 M 
and M sec = 0.15 M , respectively, with the latter value 
being an upper limit on both the companion and its 
minidisk. The fragmentation process continues even 
in the 0.7-Myr-old system and another fragment with 
mass < 0.08 M emerges in the disk of the primary at 
r ~ 200 AU. In the meantime, the outer fragment has mi- 
grated inward to approximately 720 AU. The subsequent 
evolution of the system is uncertain since we terminate 
the simulation due to enormous computational load. We 
hope that model 3 can ultimately evolve into a binary or 
triple system with low-mass ratio and (at least) one of 
the companion being in the brown dwarf mass regime. 

Three interesting features in Figure [T] are worth of 
specific emphasis. First, young protostellar disks (t < 
0.1 Myr) are conspicuously non-axisymmetric, even in 
the M c \ = 0.2 M model, but evolve with time toward 
a more regular, axisymmetric state. However, the rate 
of this evolution is much faster in models with lower 
cloud core masses. This is primarily caused by the fact 
that mass loading from the infalling envelope is the main 
driving force that sustains gravitational instability in the 
disk and the duration of this infall (i.e., the lifetime of 
the E PSF) is near-line arly proportional to the cloud core 
mass (|Vorobvovl 12010) . Second, a fairly large portion of 
the disk, irrespective of the model, is characterized by 
r > 0.1, which implies that the optical depth effect can- 
not be neglected in the early evolution of protostellar 



disks. We will return to this question in more detail in 
§ 13.21 Finally, sufficiently massive cloud cores may form 
binary or multiple systems with brown dwarf and/or sub- 
solar companions. In the future higher-resolution simu- 
lations, we hope to see planetary- mass objects forming 
in the disk. 

3.2. Disk radial structure 

Due to limited spatial resolution of modern observa- 
tional facilities at long wavelengths, vital information 
about the structure of protostellar disks is often re- 
trieved from modeling the spectral energy distribution 
based on analytical disk models that assume gas sur- 
face dens ities and temp eratu res to be power laws in ra- 
dius fe.g.lLoonev et al .112 0031 : 1 Andrews fc Williams! 120051 : 
lAndrews et al.l l2QQ9h . Our numerical simulations can 
provide a useful starting point for such an analysis by 
constructing the typical azimuthally- averaged radial pro- 
files of various physical characteristics of protostellar 
disks. 

Figure [3] presents such radial profiles for the gas sur- 
face density E (first row), midplane gas temperature T 
(second row), optical depth r (third row), enclosed mass 
M(r) (forth row), and angular velocity fl (fifth row). In 
particular, black solid lines correspond to model 1, while 
black dashed and dash-dotted lines present the data for 
models 2 and 3, respectively. The vertical columns from 
left to right correspond to four distinct times after the 
disk formation: t = 0.1 Myr, £=0.3 Myr, £=0.5 Myr, and 
t = 0.7 Myr. The vertical dotted lines show the disk 
radii r d in model 1 (left line) and model 2 (right line) 
as derived from a characteristic gas surface density of 
E d 2e = 0.1 g cm -2 for the disk to envelop e transition 
and radial gas velocity (see IVorobvov|[20Tol for details). 
The horizontal dotted lines in the third row mark charac- 
teristic optical depths of r=l (top) and r=0.1 (bottom). 
The sloped dotted lines in the bottom row show the an- 
gular velocity of a pressure-free gas in the gravitational 
field of a 5.0 M star. The red solid and dashed lines and 
also the blue dashed line are the best fits to the model 
gas density and temperature distributions. 

The top row in Figure [3] demonstrates that protostellar 
disks of equal age formed from cloud cores of greater mass 
are generally characterized by higher gas surface densities 
E, larger radii r d , and, consequently, greater masses M d . 
For instance, the 0.1-Myr-old disk in model 1 has E(r = 
20 AU) = 42 g cm" 2 , r d =117 AU, and M d = 0.03 M©, 
while the corresponding values in model 2 are E(r = 
20 AU) = 208 g cm" 2 , r d =378 AU, and M d = 0.2 M . 
This tendency is a mere consequence of an increased mass 
reservoir in the envelope, increased mass infall rate onto 
the disk, and insufficiently fast radial mass transport in 
low-mass disks dominated by viscous torques with typical 
a of the order of 0.005-0.01. As a result, low-mass disks 
cannot efficiently transport matter from the disk's outer 
edge (where most of the envelope gas lands) onto the 
forming star and the disk starts to grow both in mass 
and radius. A similar tendency of T Tauri disks being in 
general denser and hotter than those o f substellar objects 
was also found by Wi ebe et al.l (|2008[ ). 

However, the trend of more massive cloud cores to 
form denser disks becomes much less pronounced for 
cores with M c \ > 0.9 Mq. Indeed, the azimuthally- 
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Fig. 3. — Azimuthally- aver aged radial profiles of the gas surface density S (first row), gas midplane temperature T (second row), 
frequency-integrated optical depth r (third row), and angular velocity fl (fifth row). The forth row presents the enclosed mass M(r) as a 
function of radius. The vertical columns from left to right show the data at different times after the disk formation: t=0.1 Myr, £=0.3 Myr, 
t = 0.5 Myr, and t = 0.7 Myr. The solid, dashed, and dash-dotted lines present the data for the M c \ = 0.2 M© model 1, M c \ = 0.85 M© 
model 2, and M c \ = 1.7 Mq model 3, respectively. The red solid and dashed lines are the best fits to the model gas surface density and 
midplane temperature distributions in model 1 and model 2, respectively. The vertical dotted lines mark the position of the disk outer edge 
in model 1 (left lines) and model 2 (right lines). The horizontal dotted lines are the characteristic optical depths of r = 1.0 and r = 0.1. 
The sloped dotted lines in the bottom row show the angular velocity of a pressure- free gas in the gravitational field of a 5.0 M© star. 



averaged gas surface density in model 3 is only marginally 
higher than in model 2, indicating the onset of a self- 
regulatory state in which the disk grows in size and mass 
(due to continuing mass loading from the infalling enve- 
lope) rather than in gas surface density. This effect is 
caused by vigourous gravitational instability and frag- 
mentation that develop in the EPSF. Indeed, as Fig- 
ure [1] demonstrates, the M c \ — 0.2 M model 1 re- 
veals a rather weak spiral structure, while the protostel- 
lar disk in the M c \ = 1.7 M model 3 is conspicuously 
non-axisymmetric, showing a well-defined spiral struc- 
ture and multiple fragments forming in the arms. These 
fragments migrate radially inward and onto the star dur- 
ing a few tens of orbital periods (often on even shorter 
timescales), thus effectively restricting the growth of the 
disk surface density and bringing the disk back to the 
border of stability after every fragmentation episode. 

The behaviour of both the midplane gas temperature 
(second row in Figure [3]) and optical depth (third row) 
in disks of equal age is similar to that of the gas surface 
density (top row) — both T and r increase along the se- 
quence of increasing cloud core masses. As in the case 
with E, this trend appears to slow down for models with 
Af c i > 0.9 Mq. 

From the third row in Figure [3] it is seen that the in- 
ner disk regions and most fragments are optically thick 



with r > 1.0, while the outer disk regions are either 
marginally optically thin (0.1 < r < 1.0) or completely 
optically thin (r < 0.1). The fact that a substantial 
fraction of protostellar disks is optically thick may have 
important consequences for the disk mass estimates us- 
ing conventional observational techniques. Indeed, the 
forth row in Figure [3] presents the enclosed mass M(r) 
(both that of the disk and the envelope) as a function of 
radius. It is evident that a considerable fraction of the 
total mass may be concentrated in regions with r > 0.1, 
which are expected to be affected by a non-negligible op- 
tical depth. This is particularly true for the early disk 
evolution. For instance, the 0.3-Myr-old disk in model 2 
has radius r& « 615 AU and mass = 0.24 M . On 
the other hand, if we count only those regions that are 
characterized by r > 0.1 and r > 1, the corresponding 
disk masses are 0.075 M and 0.14 M , respectively, in- 
dicating that as much as 70% of the total mass content 
may be hidden from our sight. As Figure [3] demonstrates, 
both E and r express a tendency to decrease with time, 
implying that the observationally inferred disk mass are 
expected to be more accurate for disks of older age. 

The bottom row in Figure [3] reveals that the 
azimuthally- averaged angular velocity has a near- 
Keplerian profile, notwithstanding conspicuous local de- 
viations from a circular gas motion evident in Figure [TJ 
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In addition, the outer disk regions exhibit a slightly over- 
Keplerian rotation, most likely due to an additive effect 
from the gravitational field of the inner disk regions. 

It is useful to construct some functions that approxi- 
mate and summarize our model radial profiles of E and 
T for disks of distinct age. To fit our model data, we use 
the following two functions 



TABLE 2 

Best fit parameters to radial profiles of £ and T 



S = S (-^— j exp 



r + b 



(14) 
(15) 



where r is the radial distance in AU, To and Eo are 
the gas temperature (in Kelvin) and gas surface den- 
sity (in g cm -2 ) at 1 AU, respectively 6 . The form of 
Equation ([H]) is chosen to approximate E in the disk 
inner and intermediate regions by a power-law function, 
while the disk outer regions, which are usually character- 
ized by gas surface density falling off with radius faster 
than a power-law, are approximated by an exponential 
function. However, if n\ is set to zero, Equation JT3J) 
approximates the whole disk by a sole exponential func- 
tion. In the case of ni ^ 0, the parameters m and a 
express the steepness of E near the disk outer edge and 
the approximate radial position of the disk outer edge, 
respectively. In the case of m =0, the aforementioned 
parameters express the steepness and the characteristic 
radius of the radial gas surface density distribution in 
the whole disk. The form of Equation (fT5j) is motivated 
by the fact that the radial gas temperature distribution 
in the midplane of a flared disk is usually well approx- 
imated by a power-law function, while the right-hand- 
side term in Equation (fT5j) is invoked to make a smooth 
transition between the disk and the envelope, the latter 
being characterized by a background temperature Tb g . 
The parameter b represents the radial position beyond 
which the gas thermal balance is controlled by the exter- 
nal background radiation with temperature Tb g rather 
than by stellar irradiation or viscous heating. 

Red lines in the upper two rows of Figure [3] present 
our best fits to the model radial distributions of E and 
T. More specifically, the red solid and dashed lines are 
the best fits to model 1 and 2, respectively. We do not fit 
model 3, since its radial profiles are not significantly dis- 
similar to those of model 2. To fit the model data, we use 
the Marquardt-Lovenberg algorithm. Because the fitting 
functions ([T3|) and JTSJ) have too many free parameters, 
this iterative algorithm takes to many iterations and of- 
ten does not converge to reliable fits. Therefore, we vary 
manually the values of ni, 77,2, and a until the best agree- 
ment with the model data is achieved and obtain the 
best-fit values for Eo, To, m, and b. In both models, Tb g 
is set to 10 K. 

We note that Equation JT3J) is meant to approximate 
the gas surface density distribution of the disk only, while 
Equation (fT5]l does that for the midplane temperature 
of both the disk and the envelope. In addition, we fit 
the disk surface density only for r > 20 AU and ex- 
clude the inner disk regions at r < 20 AU, where a 
local peak/flattening in E is often seen. This feature 

6 The parameters a and b are usually large numbers. 
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is hardly present in the corresponding radial profiles of 
T. This local maximum or flattening in E is seen in 
many multidimensional numerical simulations of circum- 
stellar disks with accret i on on t o a forming; sta r (e.g . 
lLaughlin fc Bodenheimerl 119941 : iKrumholz et all 120071 : 
Kratter et al. 2010a). In our numerical simulations this 
feature is likely caused by an absorbing inner computa- 
tional boundary — the disk material is allowed to freely 
flow through the sink cell but is not allowed to flow out 
of it. Indeed, we ran a few models with a smaller sink 
cell (r sc =2 AU) and found that the position of the peak 
in E was always located a few AU away from the inner 
computational boundary. 

The resulting best-fit parameters to model radial pro- 
files of E and T in models 1 and 2 _are listed in Table [2j 
We find that the radial profiles of E in model 1 are best 
fitted by a sole exponential function (m = 0) rather than 
by the product of a power-law (ni 7^ 0) and an exponen- 
tial function. Indeed, the black solid lines in the top 
panels of Figure [3] show the lack of a distinct slope in 
the log E - log r space. Viscosity-dominated poly tropic 
disks are characterize d by radial gas surface de nsity pro- 
files of similar form ([Vorobvov fc B asu 2009b), suggest- 
ing that the disk in model 1 is shaped by viscous torques 
rather than by gravitational ones. The similarity solution 
for a viscous disk with v = ac s Z oc r +1 predicts the in- 
verse proportionality of E with radius ([Hartmann et al.l 
1998). The intermediate disk regions are indeed char- 
acterized by the "classic" E cx r _1 scaling, which can 
be expected from the radial temperature dependence 
T cx r -0-5 , typical for flared disks, and from the radial 
vertical scale height dependence Z cx r 125 obtained in 
detailed vertical disk structure models of D'Alessio et al. 
(1999). At the same time, the inner disk regions in 
model 1 show a shallower than E cx r _1 scaling, while the 
outer disk regions demonstrate a notably steeper scaling 
with radius. This steepening may partly be caused by 
flattening of the corresponding T profiles, making the 
sound speed c s virtually independent of radius at large 
radii, and partly by disk flaring, which in our numerical 
simulations is described by a steeper than D'Alessio et 
al.'s dependence of Z on radius at large radii (see figure 
7 in lVorobvovll2009bl ). _ 

The radial profiles of E in model 2 (and model 3) show 
a mixed behaviour. At the early evolution t < 0.5 Myr, 
they are best fitted by the product of a power-law and 
exponential function. In the subsequent evolution, how- 



Properties of embedded disks 



9 



Radial distance (AU) Radial distance (AU) Radial distance (AU) 



10 100 1000 10 100 1000 10 100 1000 



0.1 Myr " 




0.3 Myr - 


0.5 Myr " 


model 1 






















model 2 


























model 3 













10 100 1000 10 100 1000 10 100 1000 



Radial distance (AU) Radial distance (AU) Radial distance (AU) 

Fig. 4. — Azimuthally- averaged gas midplane temperature T 
(solid lines) and irradiation temperature T- irr (dashed lines) as a 
function of radius in model 1 (top row), model 2 (middle row), and 
model 3 (bottom row). The left, middle, and right columns are the 
data obtained at t = 0.1 Myr, t = 0.3 Myr, and t = 0.5 Myr after 
the disk formation, respectively. The vertical dotted lines mark 
the position of the disk outer edge in each model. 

ever, E begins to approach a pure exponential profile. 
To demonstrate this trend, we fit the radial profile of E 
at t = 0.7 Myr by both the product of a power-law and 
an exponential function (m ^ 0, red dashed line) and 
also by a pure exponential function (ni =0, blue dashed 
line). It is obvious that the latter case is more favourable. 
It is known that protostellar disks dominated by gravity 
are characteri zed by gas surface density profiles that scale 
asEocr 1 ' 5 (IVorobvov fc Basull2QQ9bl iRice at al.ll2010f ). 
while the radial profiles of E in viscosity-dominated disks 
cannot be descr ibed by a power-law func tion with a dis- 
tinct exponent ([Vorobvov fc Basu1l2QQ9b[ ). This suggests 
that the 0.5-Myr-old disk in model 2 makes a smooth 
transition from the gravity-dominated to the viscosity 
dominated stage. 

A visual inspection of Figure [3] reveals that protostellar 
disks diminish, cool, and expand radially outward with 
time. This trend is also reflected in the time behaviour 
of the fitting parameters. For instance, in both models 
Eo and To decrease with time. Concurrently, the expo- 
nents m in model 1 and n\ in model 2 decline notably 
with time, indicating a progressive shallowing of the cor- 
responding gas surface density profiles. More specifically, 
the exponent rt\ takes a value of 1.7 in the 0.1-Myr-old 
disk of model 2 (early embedded phase) and drops to 1.0 
in the 0.7-Myr-old disk (T Tauri phase). In this con- 
text, it is interesting to note that T Tauri- type disks in 
the 1.0-Myr-old Ophiuchus star forming r egions reveal E 
gradi ents with a mean value of n\ w 0.9 ([Andrews et al.l 
2009). The parameter b decreases with time, indicat- 
ing that the thermal balance in the disk outer regions 
becomes dominated by the external environment rather 
than by internal causes such as viscous (or shock) heat- 
ing or stellar irradiation. We note that while the radial 
profiles of E shallow with time, the slope of the corre- 
sponding radial profiles of T, expressed by the exponent 
n2, remains virtually independent of time. 

It is interesting to compare the azimuthally- averaged 



gas midplane temperature T with the irradiation tem- 
perature Ti rr (^-independent) in our model protostellar 
disks at different stages of the evolution. The difference 
between the two values determines the extent to which 
the disk is locally non-isothermal in the vertical direction. 
Figured] presents the radial profiles of T (solid lines) and 
Ti rr (dashed lines) in model 1 (top row), model 2 (middle 
row), and model 3 (bottom row). The vertical columns 
correspond to different disk ages: t = 0.1 Myr (left), 
t = 0.3 Myr (middle), and t=0.5 Myr (right). The verti- 
cal dotted lines mark the position of the disk outer edge 

ra- 
it is obvious that the midplane and irradiation tem- 
peratures are virtually indistinguishable in the low-M c i 
model 1, irrespective of the evolutionary stage, indicating 
that low-mass disks are nearly isothermal in the vertical 
direction. Only more massive disks formed from cloud 
cores with M c \ > 0.9 M can show a notable negative 
vertical temperature gradient in the inner disk regions 
and this phenomenon is localized to the early disk evolu- 
tion. This temperature gradient is likely caused by vis- 
cous heating, which scales as Q 2 oc r -3 , t hus operating 
preferentiall y in the disk inner regions (e.g lKratter et al.l 
120081 l2010al ). and also, to a lesser extent, by heating via 
shock waves with typical for our numerical simulations 
Mach numbers of the order of 1-2. An example of the 
latter effect is seen in model 1 near the outer edge of the 
0.1-Myr-old disk (upper- left panel), where an accretion 
shock caused by the infalling envelope creates a local 
maximum in the gas midplane temperature. We note 
that in Figured! we have averaged the gas midplane tem- 
peratures in the azimuthal direction, which may wash 
out large local azimuthal variations. We find that the 
gas midplane temperature in dense spiral arms and, espe- 
cially, in massive fragments can be considerably greater 
than Ti rr and reach several hundreds of Kelvin. 

3.3. Time evolution of disk, stellar, and envelope 
masses 

Integrated disk and envelope masses can also give us 
an insight into the rate at which protostellar disks evolve 
with time. The left column in Figure [5] presents the stel- 
lar mass M* (solid lines), disk mass (dashed lines), 
envelope mass M env (dash-dotted lines), and also the 
disk-to-star mass ratio £ (thick solid lines) as a func- 
tion of time passed since the onset of cloud core collapse. 
The right column shows the disk outer radius r^- The 
top row corresponds to model 1, while the middle and 
bottom rows show the data for model 2 and model 3, re- 
spectively. The vertical dotted lines mark the onset of the 
Class I (left) and Class II (right) phases, which are de- 
termin ed using a classification breakdown of I Andre et al.l 
(11993ft based on the mass remaining in the envelope (see 
IVorobvovll2010l for details). Note that the masses and 
disk-to-star mass ratios in the left column are plotted at 
different scales. 

The evolution of the disk-to-star mass ratio is most 
illustrative of changes occurring in protostellar disks with 
time. Soon after the disk formation in the Class or early 
Class I phase 7 , £ quickly reaches a maximum value and 

7 In fact, our disks form somewhat later than in the reality due 
to the use of the sink cell in our numerical simulations. However, 
the difference in the formation times is only of order unity. 
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Fig. 5. — Left column. Time evolution of the star mass 
(solid lines), disk mass (dashed lines), envelope mass (dash-dotted 
lines), and disk-to-star mass ratio (thick solid lines) in model 1 
(top), model 2 (middle), and model 3 (bottom). Right column. 
Time evolution of the disk radius r& (solid lines) in model 1 (top), 
model 2 (middle), and model 3 (bottom), vertical dotted lines 
mark the onset of the Class I phase (left lines) and Class II phase 
(right lines). 

then gradually declines with time. The peak disk-to- star 
mass ratio in model 1 is notably smaller (£ < 0.25) than 
in model 2 (f < 0.6) and model 3 (f < 0.7). Higher f 
favour gravitational instability and disk fragmentation, 
as reflected in disk images of Figure [TJ 

It is evident that £ declines with time much faster in 
model 1 than in model 2 and 3, reflecting a progressively 
faster loss of disk material in lower- M c \ models due to 
accretion onto the star and, to a lesser extent, due to 
viscous expansion of the disk outer regions. This effect, 
however, is a consequence of the initial conditions. The 
surface density of our initial cores (see Equation [T0|) is 
close to the integrated column density of Bonnor-Ebert 
spheres with similar form (r ou tAo = 6) and density en- 
hancement (A = 1.2), which implies that the low- mass 
cores are characterized by smaller sizes and higher den- 
sities than their high-mass counterparts (see Table [I]). 
Such low-mass cores empty their mass reservoir faster, 
as seen from the time behaviour of the envelope mass in 
Figure [5] and can be deduced from simple free-fall-time 
arguments. This, in turn, means that disks in higher- M c \ 
models are replenished with the envelope material for a 
longer time, leading to a slower decline in £. The situ- 
ation could be reversed if we had formed more massive 
cores by taking a greater positive density enhancement 
A but keeping the core size (r out ) and form fixed. In this 
case, more massive cores are expected to empty their 
mass rese rvoirs faster due to a higher infall rate onto the 
disk (e.g. IKratter et al.ll2QlQal ). 

Another interesting feature, related to the disk physics 
rather than to the initial core setup, is that the disk 



mass is steadily growing during the Class phase and, 
to a lesser extent, in the Class I phase and shows a 
trend for decline only in the Class II phase. This im- 
plies that the rate of mass transport of matter from 
the disk outer regions onto the star is smaller than 
the rate of mass loading onto the disk from the in- 
falling envelope (see also lVorobvovl l2QQ9al iBolevI 120091 : 
iKratter et al.ll2010aj ). This phenomenon appears to be a 
necessary but not sufficient condition for disk fragmen- 
tation and development o f the burst mode of accretion 
([Vorobvov fc Basul 1201 0b). Other factors such as insuf- 
ficiently large disk radii and masses may suppress disk 
fragmentation in l ow-M c i models, as found in many pre- 
vious studies (e.g.lMatzner LeviiJl2QQ5l [Rafikov 2005; 
IKratter et aTll2010b[ ). even if the rate of mass infall onto 
the disk is higher than the rate of mass accretion onto 
the star. 

The time evolution of the stellar mass in the interme- 
diate and upper- M c \ models shows episodes of (almost) 
instantaneous increase, which are a manifestation of the 
accretion burst phenomenon (note that the disk mass 
shows correlated drops). Such episodes are not present 
in the low-M c i model. The physical significance of such 
sharp increases in stellar mass for the early stellar evo- 
lution remains to be understood. 

The time evolution of the disk radius ra (right col- 
umn in Figure[5]) reveals significant differences in models 
with distinct cloud core masses. In the M c \ = 0.2 M 
model 1, the disk radius Steadily increases with time, 
reaches a maximum of rd=240 AU at t ~ 0.5 Myr and 
appears to gradually decline afterward, reflecting the on- 
going stellar accretion and disk viscous dispersal. Most 
of the disk growth occurs in the early Class II phase due 
to viscous expansion. On the other hand, the disk in 
model 2 (M c \ = 0.85 M ) and, especially, in model 3 
(M c i = 1.7 Mq) shows large scale radial pulsations, 
which are particularly strong in the Class I and early 
Class II phases. 

These episodes of disk expansion and contrac- 
tion are byproducts of the burst m ode of accretion 
(|Vorobvov fc Basul 120051 . 120061 . l20Tobh . When a frag- 
ment forms in the disk via fragmentation of dense and 
cold spiral arms, it is quickly driven onto the star via ex- 
change of angular momentum with the arms. This fast 
migration produces a transient episode of disk expansion 
due to conservation of angular momentum 8 , which is fol- 
lowed by disk contraction. Such disk pulsations pursue 
as long as the disk is capable of fragmenting and are 
usually terminated soon after the Class II phase ensues 
and the mass loading from the envelope diminishes. The 
M c \ = 0.2 model 1 is stable against fragmentation (see 
Figure [1]) and this explains why we do not see large scale 
disk pulsations in this model. 

Disk radial pulsations may be of great significance for 
the giant planet formation mechanism via direct gravi- 
tational instability. Disk contraction leads to a transient 
episode of density increase, during which the disk may 
give birth to a set of protoplanetary embryos. Most of 

8 Globally, angular momentum of the disk plus envelope system 
is not conserved since we use the sink cell through which angular 
momentum is carried away from the system by means of protostel- 
lar jets. However, mathematically, the total angular momentum 
(both t hat of the system and tha t of the sink cell) is conserved, see 
tests in I Vorobvov fe Basul ([20061 ) . 



Properties of embedded disks 



11 



them would not survive through the EPSF and are either 
driven onto the star or dispersed. However, when mass 
loading from the natal cloud core diminishes and the 
main fragmentation phase ends, the final disk contraction 
may give birth to a last and survivabl e set of gas giants 
on wid e and relatively stable orbits ([Vorobvov fc Basul 

HoToi). 

4. DISK STRUCTURE AND EVOLUTION ALONG THE LINE 
OF INCREASING CLOUD CORE ROTATION RATES 

In this section, we study the structure and time evolu- 
tion of protostellar disks formed via the gravitational col- 
lapse of cloud cores characterized by distinct rotational 
energies. For this purpose, we consider three cloud cores 
with the ratio of the rotational to gravitational energy 
f3 = 2.8 x 1(T 3 (model 4), f3 = 5.6 x 1(T 3 (model 2), 
and (3 = 2.3 x 10 -2 (model 5). The ratio /3 is varied 
by increasing the rotational energy E rot , while the other 
model parameters, in particular the cloud core masses, 
are kept identical. 

Figure [6] presents images of the gas surface density (in 
g cm -2 , log units) for model 4 (top row), model 2 (middle 
row), and model 5 (bottom row). The vertical columns 
from left to right correspond to three distinct evolution 
times since the disk formation: t = 0.1 Myr, t = 0.3 Myr, 
and t = 0.5 Myr. The meaning of arrows and yellow 
contour lines is the same as in Figure [TJ We note the 
spatial scale increases from 800 x 800 AU for models 4 
and 2 to 1800 x 1800 AU (at t = 0.1 Myr and t = 0.3 Myr) 
and 5000 x 5000 AU (at t = 0.5 Myr) for model 5. 

Figure [6] reveals similarities and also remarkable differ- 
ences in the time evolution of protostellar disks formed 
from cloud cores of distinct j3. While protostellar disks in 
all three models exhibit gravitational instability, which is 
strongest in the early evolution but notably weakens with 
time, only the f3 = 2.3 x 10 -2 model 5 shows vigorous 
disk fragmentation, with some of the fragments having 
survived to the late evolution stage. The other two mod- 
els with lower values of j3 experienced disk fragmentation 
only in the very early stage (t < 0.1 Myr, upper- left and 
upper- middle panels). However, the fragments did not 
survive, being either dispersed due to insufficient mass 
and/or numerical resolution or accreted onto the central 
star. By t = 0.5 Myr, the high-/3 model 5 appears to 
have formed a binary (or triple?) system, with masses of 
the primary and secondary being 0.45 M and 0.1 M©, 
respectively. In fact, the latter value is an upper limit 
on both the secondary companion and its disk and the 
companion itself may end up as an upper-mass brown 
dwarf. The separation between the primary and sec- 
ondary was about 400 AU in the 0.3-Myr-old disk and it 
has increased by an order of magnitude in the 0.5-Myr- 
old disk, implying that the system may eventually break 
up and give birth to a rogue upper-mass brown dwarf. 
Figures Q] and [6] demonstrate that binary /multiple sys- 
tems with low mass ratios can form as a result of disk 
fragmentation and possible progenitors are cloud cores 
of sufficiently high original mass and/or sufficiently high 
rotation energy. The latter re quirement for disk fragmen- 
tation was also reported bv iRice at al.l (j2010). High-/3 
cores seem to form binaries/multiples at a notably larger 
spatial separation. 

Our present and previous numerical simulations 
([ Vorobvov fc Basul I2006L l2010bh reveal that most frag- 



ments migrate radially inward and (possibly) onto the 
star, while only a few migrate outward. Fragments mi- 
grate preferentially inward because they form within spi- 
ral arms and loose their angular momentum via grav- 
itational interaction with part of the arm located at 
a larger radial distance. In turn, this element of the 
arm experience s a positive torqu e and moves outward 
(see figure 5 in I Vorobvov fc Basul (|2006f ) and animation 
at www.astro.uwo.ca/~vorobyov). Even if a fragment 
forms at the very end of the arm, a continuing infall of 
the envelope material exerts a negative torque on this 
fragment, driving it radially inward. This means that 
only those fragments that form near the disk outer edge 
and in the final stages of the EPSF (when infall from 
the envelope diminishes) may ultimately survive and ei- 
ther slowly migrate outward (bottom row in Figure HQ) 
or stabilize on orbits of order 100 AU ([ Vorobvov fc Basul 
[20103). 

The growing susceptibility of protostellar disks to frag- 
mentation along the line of increasing f3 can be ex- 
plained by the fact that the higher- j3 cores form disks 
of larger size and greater disk-to-star mass ratio £. 
Both effects are mostly caused by the associated in- 
crease in the centrifugal radius r c f = Q?r A /GM(r) of 
a gas layer originally located at the radial distance r 
from the center. In particular, the effect of disk radius 
on f ragmentation has been addressed by many authors 
fe.g.lRafikovll2005l:lMatzner fc Levii 3l2QQ5t IKra tter et al.l 
l2010bl : IVorobvov fc Basul l2Q10bh who have found that 
fragmentation is suppressed in small-size disks due to 
either stellar irradiation or viscous heating. 

The left column in Figure [7] presents the time evolu- 
tion of disk masses (dashed lines), stellar masses (solid 
lines), envelope masses (dash-dotted lines), and disk-to- 
star mass ratios (thick solid lines), while the right col- 
umn shows the time evolution of disk radii (solid lines). 
In particular, the top row presents the data for the 
j3 = 2.8 x 10 -3 model 4, while the middle and bottom 
rows correspond to the j3 = 5.6 x 10 -3 model 2 and 
/3 = 2.3 x 10 -2 model 5, respectively. The vertical dot- 
ted lines mark the onset of the Class I (left) and Class II 
phases (right) in each model. 

In all three models, the disk-to-star mass ratio £ 
quickly grows with time in the Class phase, attains 
a maximum value in the late Class or early Class I 
phase, and declines steadily afterwards. In the low- and 
intermediate- /? models, the disk mass is always lower 
than that of the central star, with maximum £ = 0.5 
(model 4) and £ = 0.6 (model 2). On the other hand, 
the disk mass in the high-/3 model 5 exceeds episodically 
that of the star in the Class I phase but drops quickly 
to £ < 0.4 in the subsequent evolution when the binary 
system starts to emerge. This suggests that protostel- 
lar disks with mass comparable to or greater than that 
of the host star must be statistically rare and this phe- 
nomenon quickly resolves into a binary /multiple system. 
A similar effect w as reported in numeri cal hydrodynam- 
ics simulations of IKratter et al.l (|2010af ). Moreover, such 
massive disks must be difficult to observe due to the ob- 
scuration of light by natal envelopes that are still of sub- 
stantial mass in the Class I phase. Protostellar disks 
formed from higher- f3 cores undergo radial pulsations of 
a notably greater amplitude. This increase is caused by 
a growing strength of gravitational instability and frag- 
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Fig. 6. — Gas surface density (in g cm -2 , log units) in three models with distinct ratios of the rotational to gravitational energy f3 and at 
distinct times t after the disk formation. In particular, the top row shows disk images in the f3 = 2.8 x 10 — 3 model 4, while the middle and 
bottom rows present those of the (3 = 5.6 x 10 -3 model 2 and f3 = 2.3 x 10 -2 model 5, respectively. Columns from left to right correspond 
to t = 0.1 Myr, t = 0.3 Myr, and t = 0.5 Myr, respectively. Superimposed on the disk images are the gas velocity fields (arrows) and yellow 
contour lines delineating disk regions with the frequency-integrated optical depth r > 0.1. 



mentation. In all three models, the radial pulsations ap- 
pear to subside with time. 

Finally, in Figure [8] we consider the azimuthally- 
averaged radial profiles (from top to bottom row) of 
the gas surface density, midplane gas temperature, op- 
tical depth, integrated mass, and angular velocity at 
t = 0.1 Myr (left column), t = 0.3 Myr (middle column), 
and t = 0.5 Myr (right column) after the disk forma- 
tion. In particular, the solid, dashed and dash-dotted 
lines correspond to model 4, model 2, and model 5, re- 
spectively. The meaning of the dotted lines is the same 
as in Figure [3j 

A visual inspection of Figure [8] reveals that protostellar 
disks formed from cloud cores of increasingly higher j3 are 
characterized by lower azimuthally-averaged gas surface 
density E, midplane temperature T, and optical depth 
r. At the same time, both the disk size and mass in- 



crease along the line of increasing /3. Lower E do not pre- 
vent higher- j3 disks from fragmenting, as Figure [6] demon- 
strates. On the contrary, these disks are more suscepti- 
ble to fragmentation, owing to their larger size and lower 
gas temperature. The latter effect is caused by lower 
stellar irradiation fluxes at larger radii and lower vis- 
cous heating (e.g. iMatzner fc Levin|[2005t iKratter et al.l 
l2QlQbh . A substantial fraction of protostellar disks in 
the low- and intermediate- j3 models 4 and 2 are optically 
thick, suggesting that the mass of such disks may be 
systematically underestimated using conventional obser- 
vational techniques. Sharp jumps in the integrated mass 
of model 5 at 400 AU (0.1-Myr-old disk), 1200 AU (0.3- 
Myr-old disk), and 4200 AU (0.5-Myr-old disk) reflect 
the formation and outward drift of a sub-stellar compan- 
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Fig. 7. — Left column. Time evolution of the star mass 
(solid lines), disk mass (dashed lines), envelope mass (dash-dotted 
lines), and disk-to-star mass ratio (thick solid lines) in model 4 
(top), model 2 (middle), and model 5 (bottom). Right column. 
Time evolution of the disk radius (solid lines) in model 4 (top), 
model 2 (middle), and model 5 (bottom). Vertical dotted lines 
mark the onset of the Class I phase (left lines) and Class II phase 
(right lines). 

5. CONCLUSIONS 

Using numerical hydrodynamics simulations in the 
thin-disk approximation, we model the formation and 
evolution of protostellar disks in the embedded phase 
of star formation when the disks are exposed to intense 
mass loading from natal cloud cores and are difficult to 
probe using conventional observational techniques. The 
key physical processes that are taken into consideration 
in our modeling include: disk self-gravity, viscous and 
shock heating, stellar and background irradiation, radia- 
tive cooling from the disk surface, and also self-consistent 
accretion of gas from the disk onto the forming star and 
from the infalling envelope onto the disk. We scruti- 
nize the structure and time evolution of protostellar disks 
formed from five model cloud cores with distinct initial 
masses in the 0.2-1.7 M range and ratios of the ro- 
tational to gravitational energy j3 in the 0.28-2.3 x 10 -2 
range. We find the following. 

• The time evolution of embedded protostellar disks 
proceeds from a conspicuously non-axisymmetric 
state toward a more regular, axisymmteric appear- 
ance. Simultaneously, the gas surface density re- 
duces and disks expand radially outward with time. 
However, the time scale for this transformation is 
faster in lower-mass disks formed from cloud cores 
of smaller mass. For instance, disks formed from 
cloud cores of M c \ = 0.2 M are virtually ax- 
isymmetric at t — 0.3 Myr after their formation, 
while more massive disks formed from cloud cores 
of M = 0.85 M retain a notable diffuse spiral 



structure at t = 0.5 Myr. In the non-axisymmetric 
state, massive disks are often lopsided and the gas 
flow exhibits non-circular motions, indicating local 
deviations from a circular rotation caused by spiral 
arms and forming fragments. 

• Disk fragmentation is seen in most models. How- 
ever, most of the fragments do not survive through 
the embedded phase and are either destroyed or 
driven onto the forming star. Only models with suf- 
ficiently high M c \ or ft (for instance, M c \ = 1.7 M 
and f3 = 5.6 x 10~ 3 or M d = 0.85 M and 
j3 = 2.3 x 10 ~ 2 ) reveal the formation of wide- 
separation binary /multiple systems with low mass 
ratios and (possibly) brown dwarf companions. In 
particular, systems formed from cloud cores with 
j3 > 2.3 x 10 ~ 2 may ultimately break up, deliv- 
ering freely wandering brown dwarfs to the natal 
star forming region. These rogue brown dwarfs 
may even possess some minidisks, the properties 
of which remain to be studied. 

• The fact that only sufficiently massive cores with 
sufficiently high initial rotational energy E rot form 
binary systems can in part account for a lower bi- 
nary fraction found in stars of lower mass. Low- 
mass stars tend to form from low-mass cores and 
only a small fraction of these cores with high 
enough E rot can give birth to a binary system. 
As the core mass grows (and so does the result- 
ing mass of the star), more and more cores can 
form binary systems because the requirement of 
sufficiently high initial rotational energy is relaxed 
along the line of increasing core mass. 

• Embedded protostellar disks of equal age formed 
from cloud cores of greater mass (but equal f3) are 
generally characterized by higher gas surface den- 
sities, higher midplane temperatures, larger radii, 
and greater masses. On the other hand, proto- 
stellar disks formed from cloud cores of higher f3 
(but equal M c \) are generally thinner and colder 
but larger and more massive. 

• The trend of more massive cloud cores to form 
denser and hotter disks diminishes for M c \ > 
0.9 M owing to the onset of vigourous gravita- 
tional instability and fragmentation. Gravitational 
torques efficiently transport matter in the form of 
fragments from the disk outer regions (where most 
of the envelope material lands and fragmentation 
takes place) onto the star, preventing the disk from 
growing in density and bringing it back to the bor- 
der of stability. However, the disk continues to 
grow in size (and mass) due to infall of the enve- 
lope material with high angular momentum (which 
may lead to another fragmentation episode). This 
process continues until the envelope is depleted of 
matter. 

• The positive difference between the midplane and 
irradiation temperatures is significant only in the 
inner disk regions (r < 50 AU), near the disk's 
outer edge and dense fragments, where it can reach 
a factor of two or more. The degree of vertical 
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Fig. 8. — Azimuthally- aver aged radial profiles of the gas surface density E (first row), gas midplane temperature T (second row), 
frequency-integrated optical depth r (third row), and angular velocity Q (fifth row). The forth row presents the enclosed mass M(r) as a 
function of radius. The vertical columns from left to right show the data at three different times after the disk formation: f=0.1 Myr (left), 
t=0.3 Myr (middle), and t = 0.5 Myr (right). The solid, dashed, and dash-dotted lines present the data for the ft = 2.8 x 10 -3 model 4, 
P = 5.6 x 10 -3 model 2, and (5 = 2.3 x 10 -2 model 5, respectively. The vertical dotted lines mark the position of the disk outer edge in 
model 4 (left lines) and model 2 (right lines). The meaning of other dotted lines is the same as in Figure [T] 

2.3 x 10 -2 . However, such massive disks are short- 



non-isothermality becomes progressively smaller in 
disks formed from cloud cores of smaller mass and 
as the disks evolve. 

• Protostellar disks in the embedded phase (espe- 
cially in the Class I phase) show radial pulsations, 
the amplitude of which increases along the line of 
increasing M c \ and j3. These pulsations are caused 
by gravitational instability and, to a greater extent, 
by disk fragmentation and accretion of fragments 
onto the protostar. The amplitude of radial pul- 
sations appears to diminish in the Class II phase 
when the envelope clears. 

• The disk-to-star mass ratio £ increases along the 
line of increasing M c \ and/or f3. A maximum value 
of j3 of order unity is achieved in single stars formed 
from cloud cores with high rates of rotation j3 > 



lived and by the end of the embedded phase they 
quickly evolve into a binary /multiple system. 

• A substantial fraction of a protostellar disk in the 
embedded phase may be optically thick, implying 
that observationally inferred disk masses may be 
underestimated. 

We also provide useful approximation formula to the 
radial profiles of the gas surface density and midplane 
temperature at different stages of the evolution and for 
disks of different mass. We did not see the formation 
of planetary-mass objects in protostellar disks as in nu- 
me rical simulations with a b arotropic equation of state 
by I Vorobvov fc Basul (|2010af ). However, we hope to see 
them forming in our future simulations because some of 
our low-mass fragments seem to disperse due to the lack 
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of numerical resolution at large radii where our logarith- 
mic grid diverges. Finally, it is worth noting that the 
results may be somewhat dependent on the chosen an- 
gular momentum and density profiles. We plan to vary 
the initial profiles in a future study. 
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